function [z, QzN] = TransitProb(rhoz,zbar,stdz,sIdio,nz,ns,multiple,Jen)


%To generate transition probability matrix for AR(1) with time-varying sigma
% Z(t+1) = zbar*(1-rhoz)+rho*Z(t)+sigma(t)*epsilon(t+1)
% conditional mean of Z(t+1)=zbar*(1-rho)+rho*Z(t)
% conditional volatility of Z(t+1) is sigma(t)
% sigma(t) takes two values as low and high state of uncertainty
%  

% variable definition
% z: Z(t)
% nz: number of grid points for z;
% rhoz: persistence of z
% sIdio: sigma(t) as vector of two values
% ns: number of grid points for sIdio
% stdz: the average of volatilities
% multiple: controls the boundary of z
% Jen: Jensen's effect adjustment factor, default is 2;

zmin = zbar - multiple*stdz;
zmax = zbar + multiple*stdz;
z    = linspace(zmin,zmax,nz);

% Use the method in Bloom 2009
QzN        = ones(nz,nz,ns);
% 1. for the lower bottom points of Z(t+1): P(Z(t+1)|Z(t) =
% CDF(-inf,mean(Z(1)+Z(2)|Z(t)) for given Z(t)
QzN(:,1,1) = normcdf((z(1)+z(2))/2,zbar*(1-rhoz)+rhoz*z-sIdio(1)^2/Jen,sIdio(1));
QzN(:,1,2) = normcdf((z(1)+z(2))/2,zbar*(1-rhoz)+rhoz*z-sIdio(2)^2/Jen,sIdio(2));
% 2. for the middle points of Z(t+1), P(Z(t+1)|Z(t) =
% CDF(mean(Z(i+1)+Z(1)) - mean(Z(i)+Z(i-1)|Z(t)) for given Z(t)
for j = 2:nz-1
    for q = 1:ns
        if q==1; QzN(:,j,q) = normcdf( (z(j+1)+z(j))/2,zbar*(1-rhoz)+rhoz*z-sIdio(q)^2/Jen,sIdio(q))-normcdf( (z(j)+z(j-1))/2,zbar*(1-rhoz)+rhoz*z-sIdio(q)^2/Jen,sIdio(q));end
        if q==2; QzN(:,j,q) = normcdf( (z(j+1)+z(j))/2,zbar*(1-rhoz)+rhoz*z-sIdio(q)^2/Jen,sIdio(q))-normcdf( (z(j)+z(j-1))/2,zbar*(1-rhoz)+rhoz*z-sIdio(q)^2/Jen,sIdio(q));end
    end
end
% 3. for the bottom points of Z(t+1), P(Z(t+1)|Z(t) =1-sum of probabilities
% of the rest from i to end-1;
QzN(:,end,1) = 1-sum(QzN(:,1:end-1,1),2);
QzN(:,end,2) = 1-sum(QzN(:,1:end-1,2),2);

end


